function cfdFlowCalculateGeometry3D
%--------------------------------------------------------------------------
%  Written by LiuHaHa @ NWPU, 2021 01 
%  Contact me at: liuzhikan@mail.nwpu.edu.cn
%==========================================================================
% Routine Description:
%   This function reads poly mesh files from "constant/polyMesh" directory
%   and stores them in the database ~domain
%--------------------------------------------------------------------------
global cfdFlow;

% Get info
theNumberOfElements = cfdFlow.mesh.numberOfElements;
theNumberOfFaces = cfdFlow.mesh.numberOfFaces;
theFaceNodesIndices = cfdFlow.mesh.faceNodes;
theNodeCentroids = cfdFlow.mesh.nodeCentroids;
theElementFaceIndices = cfdFlow.mesh.elementFaces;
owners = cfdFlow.mesh.owners;
neighbours = cfdFlow.mesh.neighbours;
theNumberOfInteriorFaces = cfdFlow.mesh.numberOfInteriorFaces;

% Initialize mesh member arrays
elementCentroids = cfdVectorList(theNumberOfElements);
elementVolumes   = cfdScalarList(theNumberOfElements);
faceCentroids    = cfdVectorList(theNumberOfFaces);
faceSf           = cfdVectorList(theNumberOfFaces);
unfaceSf         = cfdVectorList(theNumberOfFaces);
faceAreas        = cfdScalarList(theNumberOfFaces);
faceRf           = cfdVectorList(theNumberOfFaces);
faceRf1          = cfdVectorList(theNumberOfFaces);
faceRf2          = cfdVectorList(theNumberOfFaces);

% Process basic face geometry
for iFace=1:theNumberOfFaces    
    theNodesIndices = theFaceNodesIndices{iFace};    
    theNumberOfFaceNodes = length(theNodesIndices);
    
    % Compute a rough centre of the face    

    local_centre = cfdVector(0,0,0);
        
    for iNode=theNodesIndices
        local_centre = local_centre + theNodeCentroids(iNode,:);
    end
    local_centre = local_centre/theNumberOfFaceNodes;
    
    centroid = cfdVector(0,0,0);
    Sf = cfdVector(0,0,0);
    area = 0;
    
    % Using the centre compute the area and centoird of vitual triangles
    % based on the centre and the face nodes    
    for iTriangle=1:theNumberOfFaceNodes
        point1 = local_centre;
        point2 = theNodeCentroids(theNodesIndices(iTriangle),:);
        if(iTriangle<theNumberOfFaceNodes)
            point3 = theNodeCentroids(theNodesIndices(iTriangle+1),:);
        else
            point3 = theNodeCentroids(theNodesIndices(1),:);
        end
        local_centroid = (point1+point2+point3)/3;
        local_Sf  = 0.5*cross(point2-point1,point3-point1);
        local_area = norm(local_Sf);
        
        centroid = centroid + local_area*local_centroid;
        Sf = Sf + local_Sf;
        area = area + local_area;
    end
    centroid = centroid/area;        
    %
    faceCentroids(iFace,:) = centroid;
    faceSf(iFace,:)        = Sf;
    faceAreas(iFace)       = area;
end

% Compute volume and centroid of each element
for iElement=1:theNumberOfElements
    
    theElementFaces = theElementFaceIndices{iElement};

    % Compute a rough centre of the element
    local_centre = cfdVector(0,0,0);
    for iFace=1:length(theElementFaces)
        faceIndex = theElementFaces(iFace);
        local_centre = local_centre + faceCentroids(faceIndex,:);
    end
    local_centre = local_centre/length(theElementFaces);
    
    % Using the centre compute the area and centoird of vitual triangles
    % based on the centre and the face nodes
    %
    localVolumeCentroidSum = cfdVector(0,0,0);
    localVolumeSum = 0;
    for iFace=1:length(theElementFaces)
        faceIndex = theElementFaces(iFace);                      
        
        Cf = faceCentroids(faceIndex,:) - local_centre;
                
        faceSign = -1;
        if iElement==owners(faceIndex)
            faceSign = 1;
        end
        
        local_Sf = faceSign*faceSf(faceIndex,:);
        
        % calculate face-pyramid volume
        localVolume = dot(local_Sf,Cf)/3;
        
        % Calculate face-pyramid centre
        localCentroid = 0.75*faceCentroids(faceIndex,:) + 0.25*local_centre;
        
        %Accumulate volume-weighted face-pyramid centre
        localVolumeCentroidSum = localVolumeCentroidSum + localCentroid*localVolume;
        
        % Accumulate face-pyramid volume
        localVolumeSum = localVolumeSum + localVolume;
    end
    %
    elementCentroids(iElement,:) = localVolumeCentroidSum/localVolumeSum;
    elementVolumes(iElement)     = localVolumeSum;
end

% Process secondary Face Geometry
for iFace=1:theNumberOfInteriorFaces   
    unfaceSf(iFace,:) = cfdUnit(faceSf(iFace,:));    
    %
    own = owners(iFace);
    nei = neighbours(iFace);        
    %
    faceRf (iFace,:)    = elementCentroids(nei,:) - elementCentroids(own,:);
    faceRf1(iFace,:)    = faceCentroids(iFace,:)  - elementCentroids(own,:);
    faceRf2(iFace,:)    = faceCentroids(iFace,:)  - elementCentroids(nei,:);
end

for iBFace=theNumberOfInteriorFaces+1:theNumberOfFaces
    unfaceSf(iBFace,:) = cfdUnit(faceSf(iBFace,:));    
    %
    own = owners(iBFace);       
    %
    faceRf (iBFace,:)    = faceCentroids(iBFace,:) - elementCentroids(own,:);
    faceRf1(iBFace,:)    = faceCentroids(iBFace,:) - elementCentroids(own,:);
    faceRf2(iBFace,:)    = 0;
end

% Save and Store
meshTemp = cfdFlow.mesh;
%
meshTemp.elementCentroids = elementCentroids;
meshTemp.elementVolumes   = elementVolumes;
meshTemp.faceCentroids    = faceCentroids;
meshTemp.faceSf           = faceSf;
meshTemp.unfaceSf         = unfaceSf;
meshTemp.faceAreas        = faceAreas;
meshTemp.faceRf           = faceRf;
meshTemp.faceRf1          = faceRf1;
meshTemp.faceRf2          = faceRf2;
%
cfdFlow.mesh = meshTemp;
